Instabilities in a free granular fluid described by the Enskog equation 
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A linear stability analysis of the hydro-dynamic equations with respect to the homogeneous cooling 
state is carried out to identify the conditions for stability as functions of the wave vector, the 
dissipation, and the density. In contrast to previous studies, this description is based on the results 
derived from the Enskog equation for inelastic hard spheres [V. Garzo and J. W. Dufty, Phys. 
Rev. E 59, 5895 (1999)], which takes into account the dependence of the transport coefficients on 
dissipation. As expected, linear stability shows two transversal (shear) modes and a longitudinal 
, ("heat") mode to be unstable with respect to long enough wavelength excitations. Comparison 

. with previous results (which neglect the influence of dissipation on transport) shows quantitative 

' discrepancies for strong dissipation. 



I. INTRODUCTION 



A simple way of capturing the dynamics of granular media under rapid flow conditions is through an idealized fluid 
of smooth, inelastic hard spheres. The inelasticity of collisions is only accounted for by a (constant) coefficient of 
normal restitution < a < 1 that only affects the translational degrees of freedom of grains. Despite the simplicity of 
the model, it has been widely used as a prototype to understand some of the physical mechanisms involved in granular 
flows, especially those related to the inelasticity of granular collisions. In particular, one of the most characteristic 
Ch features of granular fluids, as compared with normal fluids, is the spontaneous formation of velocity vortices and density 
' clusters when evolving freely. This clustering instability can be well described through a linear stability analysis of 
the hydrodynamic equations and follows from the presence of a dissipation term in the equation for the balance of 
energy. An important feature of this instability is that it is confined to long wavelengths (small wave numbers) and 
O | so it can be avoided for small enough systems. First detected by Goldhirsch and Zanetti [lj and McNamara Q, the 
clustering problem has attracted much attention in the past few years, especially from a computational point of view 

^ ■ In the case of a low-density gas, accurate predictions for the unstable hydrodynamic modes have been made from 
r — ' the (inelastic) Boltzmann equation 0,0,0- In particular, a critical length L c is identified, so that the system becomes 
On : unstable when its linear size is larger than L c . The dependence of L c on the coefficient of restitution a predicted 
by kinetic theory compares quite well with numerical results obtained by using the direct simulation Monte Carlo 



method 0. 

ly-j ■ For finite higher densities, these instabilities have been studied by several authors using macroscopic or kinetic 
equations A careful study of the dispersion relations has been recently carried out by van Noije and Ernst 

1*-; ■ 9] from the Enskog kinetic theory but neglecting any dependence of the pressure and of the transport coefficients on 
inelasticity. Specifically, they assume that the expressions for the hydrostatic pressure, the shear viscosity, the bulk 
viscosity, and the thermal conductivity are the same as those for the elastic gas |lfj) , except that the temperature still 
depends explicitly on time to account for the homogenous cooling. However, given that the effect of inelasticity on 
dense fluid transport is quite important in the undriven case [Tllll2 | (see for instance, Fig. ^ below) the assumptions 
made in Ref. 9] could be only justified for very small dissipation (say for instance, a 0.99). 

Although the predictions made by van Noije and Ernst [j| compares reasonably well with molecular dynamics 
simulation results, it is worth to assess to what extent the previous results are indicative of what happens when 
the improved expressions for the inelastic transport coefficients are considered For this reason, in this paper 

1 revisit the (unstable) hydrodynamic-mode problem of a granular fluid described by the inelastic Enskog kinetic 
theory ^| . In spite of the explicit knowledge of the Enskog transport coefficients ^l] , I am not aware of any previous 
solution of the linearized hydrodynamic equations for a moderately dense granular gas. 

The Enskog kinetic equation can be considered as the extension of the Boltzmann equation to finite densities. 
As happens for elastic collisions, the inelastic Enskog equation provides a semiquantitative description of the hard 
sphere system that neglects the effects of correlations between the velocities of the two particles that are about to 
collide (molecular chaos assumption). The Enskog approximation is expected to be valid for short times since as the 
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system evolves corrections to the Enskog equation due to multiparticle collisions, including recollision events ("ring" 
collisions) should be incorporated. The latter are expected to be stronger for fluids with inelastic collisions where 
the colliding pairs tend to be more focused. Therefore, some deviations from molecular chaos have been observed 
in molecular dynamics (MD) simulations [l4| of granular fluids as the density increases. Although the existence 
of these correlations restricts the range of validity of the Enskog equation, the latter can be still considered as a 
good approximation especially at the level of macroscopic properties (such as transport coefficients). In particular, 
the Enskog results presents quite a good agreement with MD simulations and even with real experiments. In the 
case of computer simulations, comparison between the Enskog theory and MD simulations in the case of the self- 
diffusion coefficient |T^ | and kinetic temperatures in a granular mixture |l6j | have shown good agreement for all a at 
n* = na 3 < 0.25 and for all densities at a > 0.9. Here, n is the number density and a is the diameter of spheres. 
More recen t ag reement has been found in the case of granular mixtures under shear flow |l7| . The Enskog transport 
coefficients [ll| have also been tested against real NMR experiments of a system of mustard seeds vibrated vertically 
[Lsl. fli^ . The averaged value of the coefficient of restitution of the grains used in this experiment is a — 0.87, which 
lies outside of the quasielastic limit (a s» 0.99). Comparison between theory and experiments (see for instance, Figs. 
10-13 of Ref. 0|) shows that the Enskog kinetic theory ^l| successfully models the density and granular temperature 
profiles away from the vibrating container bottom and quantitatively explains the temperature inversion observed in 
experiments. All these results clearly show the applicability of the Enskog theory for densities outside the Boltzmann 
limit (n* — ► 0) and values of dissipation beyond the quasielastic limit. In this context, one can conclude that the 
Enskog equation provides a unique basis for the description of dynamics across a wide range of densities, length scales, 
and degrees of dissipation. No other theory with such generality exists. 

The explicit knowledge of the Navier-Stokes transport coefficients as well as of the cooling rate for inelastic hard 
spheres [Tl| allows one to solve the linearized hydrodynamic equations around the homogeneous cooling state (HCS) 
and identify the conditions for stability as functions of the wave vector, the dissipation, and the density. In the 
low-density limit, previous results derived from the Boltzmann equation are recovered Q. Linear stability analysis 
shows two transversal (shear) modes and a longitudinal ( "heat" ) mode to be unstable with respect to long wavelength 
excitations. The corresponding critical values for the shear and heat modes are also determined, showing that the 
clustering instability is mainly driven by the transversal shear mode, except for quite large dissipation. As expected, 
these results agree qualitatively well with those previously derived in Ref. 0. On the other hand, at a quantitative 
level, the comparison carried out here shows significant differences between both descriptions as the collisions become 
more inelastic. 

The plan of the paper is as follows. In Sec. [HJ the basis of the hydrodynamic equations for a dense gas as derived 
from the (inelastic) Enskog equation is described. The explicit dependence of the transport coefficients and the cooling 
rate on dissipation is also illustrated for some values of density to show that the influence of inelasticity on transport 
is in general quite significant. Section IlIII is devoted to the linear stability analysis around the HCS and presents the 
main results of this paper. The paper is closed in Sec. IIVI with some concluding remarks. 



II. HYDRODYNAMIC DESCRIPTION 

We consider a granular fluid composed of smooth inelastic hard spheres of mass m and diameter a. Collisions 
are characterized by a (constant) coefficient of normal restitution < a < 1. At a kinetic level, all the relevant 
information on the system is given through the one-particle velocity distribution function, which is assumed to obey 
the (inelastic) Enskog equation |13| . From it one can easily get the (macroscopic) hydrodynamic equations for the 
number density n(r,t), the flow velocity u(r,t), and the local temperature T(r,t) 

D t n + nV- u = , (1) 



pAu + VP = 0, (2) 



£> t T+|-(V-q+P:Vu) = -CT. (3) 

In the above equations, D t = dt + u • V is the material derivative, p — nm is the mass density, P is the pressure 
tensor, q is the heat flux, and C is the cooling rate due to the energy dissipated in collisions. The practical usefulness 
of the balance equations CEJ-© is limited unless the fluxes and the cooling rate are further specified in terms of the 
hydrodynamic fields and their gradients. The detailed form of the constitutive equations and the transport coefficients 
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appearing in them have been obtained by applying the Chapman-Ensk og m ethod [Toj to the Enskog equation. To 
first order in the gradients, the corresponding constitutive equations are [ljj 



Pij = pSij - rj ( VjUi + ViUj - -Sij V • u ) - 7<5y V • u, 



(4) 



-kVT — /uVn, 



(5) 



C = Co + CiV 



(G) 



Here, p is the hydrostatic pressure, 77 is the shear viscosity, 7 is the bulk viscosity, n is the thermal conductivity, 
and fi is a new transport coefficient not present in the elastic case. The expressions for the pressure, the transport 
coefficients and the cooling rate can be written in the forms 



p = nTp*(a, 



k = K Q n*(a, 0), // = 4>), Co = —Co («i 

n 770 



(8) 



where 7/0 = 5(?7jT) 1 / 2 /16cr 2 7r 1 / 2 and ko = 157/o/4777 are the low-density values of the shear viscosity and the thermal 
conductivity in the elastic limit, respectively. The quantities p* , 77*, 7*, ft*, /j*, Q, and Ci are dimcnsionless functions 
of the coefficient of restitution a and the solid volume fraction </> = irna 3 /6. Their explicit expressions are given in 
Appendix [3] and more details can be found in Ref. For elastic collisions (a — 1), 0), Co (1, </>)*, an d Ci(l) 4>) 
vanish, while the expressions of 77* (1, <j>), 7*(1, <fi), and k*(1, <j>) coincide with those obtained for a dense gas of elastic 
hard spheres ^(j- I n the low-density limit (0 = 0), 7* (a, 0) = (i(a, 0) = and the results derived for a dilute granular 
gas are recovered 0. As pointed out before, the new transport coefficient /1 is not present for elastic collisions and 
may play an important role to accurately describe some situations of real granular materials, such as a "temperature 
inversion" observed in vibrofludized systems [19(. 

The reduced quantities r)*(a,4>)/r)*(l,(f)), n*(a, 4>)/n*(l, (f), /j*(a,(/>), and (i(a,4>) are plotted in Fig. functions 
of the coefficient of restitution for three different values of the solid volume fraction (p. As said in the Introduction, 
these quantities, along with the pressure, were assumed to be the same as in the elastic case in the stability analysis 
carried out in Ref. 0, i.e., p*(a, — * p*(l, </>), 77*(a, 4>) — > 77* (1, </>), n*(a, — > 0), and H*(a, </>) = Ci( a 7 0) — * 0. 
Figure ^ shows that in general the influence of dissipation on the transport coefficients and the cooling rate is quite 
significant and so their functional form differs appreciably from their elastic form. This means that the predictions 
made in Ref. might quantitatively differ from those obtained here as the rate of dissipation increases. This will 
be confirmed later. We also see that, for a given value of a, rf{a.,4>)/rf{l,4>) and K*{a,4>)/ K*{\,(j>) decrease as the 
density increases, while the opposite happens in the cases of p,*(a,(f>) and |Ci( a ) 4>)\- 

When the expressions of the pressure tensor, the heat flux and the cooling rate are substituted into the balance 
equations Q~0 one gets the corresponding Navier-Stokes (closed) hydro dynamic equations for n, u and T. They 
are given by 



D t n + nV ■ u = 0, 



(9) 



D t Ui + (77,777) v%p = (77,771) Vj 



77 ViUj + VjUi - -SijV 



-7%V 



(10) 



(A + Co)T + ^- P X7 ■ u = ^-V • (kVT + /iVn) 

on 077 

2 

3n 



77 ( ViUj + VjUi - -5ijV • u ) + 7<%V • u 



- TCiV • u. 



(11) 



Note that consistency would require to consider up to second order in the gradients in the expression (JSJ for the cooling 
rate, since this is the order of the terms in Eq. coming from the pressure tensor and the heat flux. However, it has 
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FIG. 1: Plot of r,*{aA)h*{l,<i>), «*(<*, *)/«*(!, <P), 0), 
different values of the solid volume fraction cj>: (a) <j> = 0, (b) (f> 



and £i(a,(f>) versus the coefficient of restitution a for three 
= O.f , and (c) <j> = 0.2. Note that & = at <j> = 0. 



been shown for a dilute gas that the contributions from the cooling rate of second order are negligible as compared 
with the corresponding contributions from Eqs. 4J. It is assumed here that the same holds in the dense case. 

The form of the Navier-Stokes equations |§)l- HII|) is the same as for a normal fluid, except for the presence of the 
contributions to the cooling rate Co and Ci an d the new transport coefficient /i in the energy balance equation. Of 
course, as Fig. ^clearly illustrates, the values of the transport coefficients are quite different, depending on the value 
of the coefficient of restitution a. 



III. LINEAR STABILITY ANALYSIS 



The hydrodynamic equations JIJ-© admit a simple solution which corresponds to the so-called homogeneous 
cooling state (HCS). ft describes a uniform state with vanishing flow held and a temperature decreasing monotonically 
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in time, i.e., 

T(t) = ™ . . (12) 

(l + Co(0)t/2) 2 

Nevertheless, computer simulations P, |2(], have shown that the HCS is unstable with respect to long enough 
wavelength perturbations. To analyze this problem, it is adequate to perform a stability analysis of the nonlinear 
hydrodynamic equations (|§|)- (|11[) with respect to the homogeneous state for small initial excitations. The lineariza- 
tion of Eqs. l|§ll- (|ll|) about the homogenous solution yields partial differential equations with coefficients that are 
independent of space but depend on time since the reference (homogeneous) state is cooling. This time dependence 
can be eliminated through a change in the time and space variables and a scaling of the hydrodynamic fields. 

We assume that the deviations Sy a (r,t) = y a (r,t) — yna(t) are small, where, 8y a (r,t) denotes the deviation of 
{n, u, T, } from their values in the homogeneous state, the latter being denoted by the subscript H. The quantities 
in the HCS verify 

X7n H =X7T H = 0, u„=0, d t lnT H = -Cwr- (13) 

To recover the results found in the dilute gas case when <j> — ► 0, 1 consider the same time and space variables as those 
used in £|, namely, 

r=\f Mt'W, i = - 2 ^r, (14) 

where Vu(t) = ^§- n H cP'i^^vu (i) is an effective collision frequency and vn(t) = y/Tjj{i)/m. Note that Vn(t) = 
tihTh /t]h0-, 0) is an effective collision frequency associated with the elastic (a = 1) shear viscosity of a dilute gas 
(0 = 0). The dimensionless time scale r is the time integral of the average collision frequency and thus is a measure 
of the average number of collisions per particle in the time interval between and t. The unit length Ujj(t)/i/jj(t) 
introduced in the second equality of (|14fl is proportional to the time-independent mean free path of gas particles. 
A set of Fourier transformed dimensionless variables are then introduced by 

Pk(r) = , w k (r) = — , 6» k (r) = , (15) 

n H v h (t) T h (t) 

where 5y^ a = {Sn-^, Wk(r), O^t)} is defined as 

5y ka (r)= J d£e-* k£ 6y a (£,T). (16) 

Note that in Eq. ||TSJ| the wave vector k is dimensionless. In terms of the above variables, the transverse velocity 
components = Wk — (wk • k)k (orthogonal to the wave vector k) decouple from the other three modes and hence 
can be obtained more easily. Their evolution equation is 

#- - Co + W k ") w k± = 0, (17) 



dr ^ ' 2 

where it is understood that ^o an( i V* are evaluated in the HCS. The solution to Eq. I|17|) is 

w k ±(k,r) = w k j_(0)exp[sj_(fc)T], (18) 

where 

s±(k)=C -lv*k 2 . (19) 

This identifies two shear (transversal) modes analogous to the elastic ones. According to Eq. Ijl9|l . there exists a 
critical wave number k s given by 

z so 



ks = y-f) ■ (20) 

This critical value separates two regimes: shear modes with k > k s always decay while those with k < k s grow 
exponentially. 
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FIG. 2: Dispersion relations for a granular fluid with a = 0.8 and <f> = 0.2. From top to bottom the curves correspond to the 
two degenerate shear (transversal) modes and the remaining three longitudinal modes. The dashed lines correspond to the 
results obtained for the shear (s±) and heat (S3) modes from the approximations made in Ref. [IJ. Only the real parts of the 
eigenvalues is plotted. 




The remaining (longitudinal) modes correspond to p^, 8^, and the longitudinal velocity component of the velocity 
field, w k M = wit • k (parallel to k). These modes are coupled and obey the equation 



dr 



(21) 



where 5y^ a (T) denotes now the set {/?k, } and M is the square matrix 



M 



-ik 

-2Qg-^*P -Q-l K *k 2 + |Ci) 

-ikp*G p -ikp* Co " ^v*k 2 ~^l*k 2 



(22) 



As before, it is understood that p*, rj*, 7*, n* , pf, Co> and Ci are evaluated in the HCS. In addition, the quantities 
g(4>) and C p (a, 4>) are given by 



d 

g{4>) = l + <j,— \n X {4>) 



(23) 



C p {a,<f>) = l+g(0) 



P*(a,^)-i 
p*{a,cj)) 



= 1 + 9(4)- 



9{<t>) 



(24) 



l + 2(l + a)teW 

where in the last equality use has been made of the explicit expression of p* given by Eq. (|A1I) . In Eqs. ill'- ill and i|24[l . 
x(</>) is the pair correlation function at contact. In kinetic theory calculations, the value of \ f° r the pre-collisional 
distribution is used, which is well approximated by local equilibrium. There are nonequilibrium corrections that can 
be calculated from the ring collision operator [2i|. However , th ese corrections are very hard to calculate and so for 
simplicity I take here the Carnahan-Starling approximation |23| : 



x{4>) 



2(1-0)3' 



(25) 



In the limit <j> — > 0, p* = g = C p = 1, 7* = £i = and Eqs. ljl"7|) - (|2*2l reduce to those previously derived for a dilute 
gas 0. 

The longitudinal three modes have the form exp[s„(fc)r] for n — 1,2,3, where s n (k) are the eigenvalues of the 
matrix M, namely, they are the solutions of the cubic equation 
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5 7 



5 , 

o7 



V C p + -77 Co + 7p Co 



~P*(2p*+3Ci)-^«*Co 



-(K*C p -p*)k 2 + Co(C P -2g) 



k 2 = 0. 



(26) 
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FIG. 3: Plot of the the critical wave numbers k s (a, <f>) and fch(a, <f)) as functions of the coefficient of restitution a for two values 
of the solid volume fraction <j>: (a) cj> = and (b) 4> = 0.2. 



As happens for a dilute gas jj] , for given values of a and (j> the analysis of Eq. l|26fl shows that at very small k all modes 
are real, while at larger k two modes become a complex conjugate pair of propagating modes. Thus, the physical 
meaning of the longitudinal modes is different from that in the elastic fluids, even when a — > 1 
The solution to Eq. (|26l) can be obtained for small k by a perturbation expansion as 

S „(fc) = ,s(°)+fc,s( 1) +fc 2 .€ ) + -- - ■ (27) 
Substituting this expansion into Eq. (|27>|) yields sf^ = 0, = —Co, s^ ' = Co, 

S « = 0, (28) 



(2) _ P 
'1 — /■* 
SO 



(C P - 2g) , 



(29) 



,(2) _ 



,2p* + 3(Ci + 2<7) 5 



6Co 



— K 

4 



(30) 



(2) ,2p*+3(Ci -2g) + 6C p 2 . 1 

n ' 

6C * 3 2 



4 = -P —* - o 7 ? - o7 • ( 31 ) 



In the case of a dilute gas (cf> — > 0), the eigenvalues s n (fc) behave as 

1 



si(fc)->-^fc 2 , (32) 



« 2 (fc) -> -Co + =>r 



3Co* ~ 4 K 



fc 2 , 



(33) 



*3(*) - Co - UtT 



1 

3Co* 



3* 



(34) 
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FIG. 4: Plot of the ratios k 3 /ko s and kh/koh as functions of the coefficient of restitution a for three values of the solid volume 
fraction cj>: (a) = 0, (b) 4> — 0.1 and (c) (j> — 0.2. Here, fcos and koh are the critical wave numbers obtained from the 
approximations made in Ref. [jj. 

Since the Navier-Stokes hydrodynamic equations are valid to second order in k, the solutions 1371) (|31|l are relevant 
to the same order. 

As said in the Introduction, although the limitations of the Enskog theory are greater than for elastic systems, 
comparison with MD simulations [TH Il6l] indicate that it is still accurate for <j> up to about 0.15 and for a greater 
than about 0.5. For higher densities the a range is more limited, but even then it captures the relevant qualitative 
features. For this reason, to illustrate the influence of both density and dissipation on instabilities, densities in the 
interval < (j> < 0.2 for 0.5 < a < 1 will be considered. 

The dispersion relations s n (k) for a fluid with a — 0.8 and <j> — 0.2, as obtained from Eq. (|20|l and the solutions of 
the cubic equation l|26|) . are plotted in Fig. [3 Only the real part (propagating modes) of the solutions to Eq. 126|) is 
represented. For comparison, the results derived for the shear (s±) and the longitudinal "heat" (S3) modes from the 
approximations made by van Noije and Ernst are also plotted. These curves can be formally obtained from the 
results derived in this paper when one takes fi* = £i = 0, and p* , 77*, 7*, and k* are replaced by their values in the 
elastic limit [Eqs. I|A18|) — (|A20|) ] . We observe that the agreement between both sets of results is in general good, the 
heat mode showing more quantitative discrepancies. Figure [3 also shows that the heat mode is unstable for k < kh, 
where kh can be obtained from Eq. (|26[1 when s = 0. The result is 



The dependence of the critical values k s and kh on dissipation is illustrated in Fig [3] for two values of cj). For a 
given value of the coefficient of restitution a, in general the corresponding critical values decrease with increasing 
density. However, there is a small region of values of a > 0.82 where the opposite happens in the case of k s . All 
the above trends are also captured by the results obtained in Ref. |9(, although quantitative discrepancies between 
both descriptions appear as the dissipation increases. To illustrate such differences, the ratios k s /k 0s and kh/k h are 
plotted versus a in Fig. 0] for different values of <f>. Here, k$ s and fcoft, are the critical wave numbers obtained from the 
approximations made in Ref. |9j. Significant differences between both analyses are clearly shown in Fig.^ especially 
for strong dissipation and moderate densities. Thus, for instance, for </> = 0.2 and a = 0.8 the discrepancies between 
both approaches for k s and kh are about 5% and 17%, respectively, while for <j> = 0.2 and a = 0.5 the discrepancies 
are about 12% and 48%, respectively. A more significant qualitative difference between our results and those obtained 
111 Ref. @ appears in the case of the ratio k s /kh- This quantity measures the separation between both critical modes. 
According to Eqs. (|20() and l|35(l . this ratio is independent of a when one neglects the influence of dissipation on the 




(35) 
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FIG. 5: Ratio k s /kh versus the coefficient of restitution a for three values of the solid volume fraction <f>: (a) (j> — 0, (b) <f> = 0.1, 
and (c) <f> = 0.2. 




FIG. 6: Dependence of ao on the solid volume fraction <j>. Points above (below) the curve correspond to systems where the 
instability is driven by the shear (heat) mode. 



pressure and the transport coefficients. However, the present results predict a complex dependence of k s /kh on a. To 
illustrate it, the ratio k s /kh is plotted versus the coefficient of restitution a in Fig. [5] for different values of density. 
It is apparent that in general both critical values k s and kh are well separated, especially for small inelasticity. The 
results also show that the influence of dissipation on the ratio k s /kh is less significant as the system becomes denser. 
In addition, for a given value of cj>, there exists a value of the coefficient of restitution (Xq(4>) for which kh > k s for 
values of a < ao- The dependence of ao on the solid volume fraction is plotted in Fig. El It is apparent that the 
value of ao decreases with increasing density. However, given that the values of ao are quite small, one can conclude 
that in practice (a > 0.375) the instability of the system is driven by the transversal shear mode since k s > kh for 
a > a (4>). 

According to these results, for not quite extreme values of dissipation (a > 0.375), three different regions can be 
identified. For k > k s all modes are negative and the system is linearly stable with respect to initial perturbations 
with wave number in this range (short wavelength region). For kh < k < k s , the shear mode is unstable while the 
heat mode is linearly stable. In this range the density (coupled to the heat mode) is also stable and so, density 
inhomogeneities can only be created due to the nonlinear coupling with the unstable shear mode pl|. Finally, if 
k < k s first vortices and then clusters are developed and the final state of the system is strongly inhomogeneous. A 
more detailed analysis of the evolution of the granular gas can be found in Ref . |3j . 

In a system with periodic boundary conditions, the smallest allowed wave number is 2it/L, where L is the largest 
system length. Hence, for given values of inelasticity and density, we can identify a critical length L c so that the 
system becomes unstable when L > L c . The value of L c is determined by equating 

^ = max{k s ,k h }, L* = t^-L c . (36) 

In Fig. [7| we show L c /\q as a function of a for different values of the solid volume fraction <fi. Here, Ao = 
[y/^^na 1 x) 1 — (5V27TX/16) 1 vh/vh is the mean free path of a hard-sphere dense gas. In all of these systems, 
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FIG. 7: The critical size L c in units of the mean free path Ao as a function of the coefficient of restitution a for three different 
values of the solid volume fraction <fr. (a) = 0, (b) <j> — 0.1, and (c) (f> — 0.2. In each case, the system is linearly stable for 
points below the corresponding curve. 



k s > kh and so 



For a given value of a, we see that the critical size (in units of the mean free path) increases with density. As a 
consequence, larger systems are required to observe the shearing instability as the fluid becomes denser. 



A well-known feature of rapid granular flows is the instability of the homogeneous cooling state against long 
wavelength spatial perturbations, leading to cluster and vortex formation. Although the origin of this instability has 
been widely explored by using computational tools, prior analytical work on this subject has been limited to weak 
inelasticity or very dilute regime. In order to gain some insight into the influence of both density and dissipation 
on the stability of the HCS, a kinetic theory description has been adopted. For moderate densities, the inelastic 
Enskog equation ^| can be considered as a valuable tool for granular media. As in the case of elastic collisions, 
the Enskog equation takes into account spatial correlations through the pair correlation function but neglects the 
velocity correlations between the particles that are about to collide (molecular chaos). The latter assumption has 
been clearly shown to fail for inelastic collisions as the density increases , so that the limitations of the Enskog 
description are greater than for elastic collisions. Due to this molecular chaos breakdown, some authors conclude 
that the Enskog equation can be insufficient to compute average properties of inelastic fluids, except for very weak 
dissipation. Nevertheless, this conclusion contrasts with previous comparisons made with MD simulations [l5lllra ll7| 
and with real experiments |18| where, at least for the problems studied there, velocity correlations do not seem to play 
an important role and the Enskog equation provides quite good estimates for the transport properties of the system. 
It is remarkable that its accuracy is not restricted to the quasielastic limit since it covers values of moderate density 
(0 < 4> < 0.15) and lar ge v alues of dissipation (0 < a < 0.5). It is possible that for situations more complex than 
those analyzed in Refs. |l5l llii ITtI fTsl. USfy . velocity correlations become important and the Enskog theory does not 
give reliable predictions. In this case, new kinetic theories incorporating the effect of velocity correlations are needed 
to describe granular flows. However, so far there is no alternative to the Enskog theory for finite density systems at 
this point. Hence, it is the most accurate theory to describe systems of interest in simulations and experiments. 

In this paper I have used the inelastic Enskog kinetic theory to perform a linear stability analysis of the hydrody- 
namic equations and identify the conditions for stability in terms of dissipation and density. The analysis is based on 
a previous derivation of the expressions of the Enskog transport coefficients and the cooling rate that, a priori, 
is not limited to small dissipation. This is the main new ingredient of this work since previous studies 9] on linear 
stability analysis for dense granular gases considered weakly inelastic systems and so, thermodynamic and transport 
properties were assumed to be the same as those of elastic hard sphere fluids. However, this assumption is expected 
to fail as dissipation increases since the form of the inelastic transport coefficients clearly differs from their elastic 
counterparts, as shown for instance in Fig. ^ 




(37) 
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The study reported here extends to higher densities a previous linear stability analysis performed for a dilute gas 
. In general, the findings agree qualitatively well with previous results [t| , showing that the effect of dissipation on 
transport coefficients do not significantly modify the qualitative form of the dispersion relations. Specifically, linear 
stability analysis shows two unstable modes: a transversal shear mode and a longitudinal "heat" mode. The instability 
of both modes is a long wavelength instability. The analysis of the dependence of the corresponding critical shear 
wave number k s [defined in Eq. H2()|l ] and heat wave number kh [defined in Eq. Ij35|l ] shows that, except for extreme 
values of dissipation, the instability is driven by the transversal shear mode. The range of values of the coefficient of 
restitution a for which kh > k s is shortened as the gas becomes denser. Thus, for </> > 0.389, k s > kh for any value of 
a. 

On the other hand, as expected, quantitative discrepancies between our results and those given in Ref. become 
significant as the dissipation increases. In particular, at a given value of density, the critical wave numbers k s and 
kh are in general underestimated (except in the case of k s for a low-density gas) when one neglects the influence 
of inelasticity on transport Q while the ratio k s /kh (which is independent of the coefficient of restitution a in 9]) 
presents a complex dependence on the rate of dissipation, as is illustrated in Fig.0 Therefore, although the description 
made by van Noije and Ernst 0] predicts reasonably well the dispersion relations as well as the long-range structure, 
one expects that the results reported here improve such predictions when one considers values of the coefficient of 
restitution a for which transport properties are clearly affected by the rate of dissipation. 

As said in the Introduction, in the case of a dilute gas (</> = 0) comparison with direct Monte Carlo simulation of 
the Boltzmann equation has shown the accuracy of the stability analysis performed in Ref. 4] . Given that the results 
reported here extends the above description to high densities, comparison with MD simulations becomes practical. In 
this context, it is hoped that the description reported here stimulates the performance of such computer simulations 
to characterize the onset and evolution of the clustering instability. As in the Boltzmann case 0, 01 > one expects 
that the Enskog results describes accurately the first stages of evolution. 

Finally, it must noted that all the results obtained in this paper has been made in the context of a very simple 
collision model where the coefficient of restitution is constant. Recent results [2^ derived with an impact- velocity- 
dependent coefficient of restitution shows that structure formation occurs in free granular gases only as a transient 
phenomenon, whose duration increases with the system size. 
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cooling rate to first order in the gradients. Partial support of the Ministerio de Educacion y Ciencia (Spain) through 
Grant No. FIS2004-01399 (partially financed by FEDER funds) is acknowledged. 



In this Appendix, the expressions for the hydrostatic pressure, the transport coefficients and the cooling rate used 
in Eqs. are given. The reduced hydrostatic pressure p* is 



The reduced transport coefficients rj*, 7*, k* , and /1* defined through the relations Q) and © are given, respectively 
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APPENDIX A 



p* =l + 2(l + a)d> X 



(Al) 



by 



V* =V* k l + -0x(l + a) +-7* 



(A2) 




(A3) 




(A4) 
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Here, the superscript k denotes the contributions to the transport coefficients coming from the kinetic parts of the 
fluxes UJ1. These kinetic contributions are 



Vk =\v„- 2<o 



l--(l + a)(l-3a)fo 
5 



(A6) 



x 1 1 + [1 + (1 + a)(t> X ] c + |fo(l + a) 2 



2a - 1 



1 + a 



3(1 + a) 



(A7) 



/4 = 2(2<-3C *)- 1 {(l + ^ln X )Co< 

+ip*(l + 0S lnp*)c- -0x(l+a) 
3 5 

x fl + ^0<V n x) [a(l - a) 
+ \ (| + a(l-a))c 
In these expressions, c, £q, i/*, and u* are functions of a and given by 



(A8) 



32(l-q)(l-2q 2 ) 
81 - 17a + 30a 2 (l - a)' 



(A9) 



(A10) 



V n = X 



I- -(I- a f 



V 64/ 



(All) 



, 33., , 19 -3a 

1 H (1 — a) H c 

16 V ' 1024 



(A12) 



Furthermore, in three dimensions the Carnahan-Starling approximation [23J for the the pair correlation function at 
contact x(0) is given by 



2(1-0)3- 

The coefficient £i appearing in the expression © for the cooling rate £ is 



(A13) 



<i 



— ( 1 + — c I Cf - 2 
32 I 64 ' C 



^x(i-« 2 ), 



(A14) 



where 126 



±A + (l + a)(i-a) 



, y *_| ( l_ a2 )(l + ^ c ) + |c( 1 + ^ c ) (1 _ a2) - 



(A15) 



241 - 177a + 30a 2 - 30a 3 + ^ (30a 3 - 30a 2 + 2001a - 1873) 



(A16) 
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3 r c 

A = -(1 + a) (1 - a)(5a 2 + 4a - 1) + — (159a + 3a 2 - 19 - 15a 3 ) 



(A17) 



For elastic collisions (a = 1), c = Q = [i* = Ci =0 and 



p*(l,0) = 1 + 4*6 7*(1,0) 



256 ,o 

OTT 



(A18) 



( 




2 



X 



-1 



(A19) 



( 




2 



(A20) 



These expressions coincide with well-know results derived for normal hard-sphere fluids [l 
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